This PR implements a revised and ―possibly― improved contact model originally proposed in the following manuscript:

Exponential Integration for Efficient and Accurate Multi-Body Simulation with Stiff Viscoelastic Contacts
Hammoud, Olivieri, Righetti, Carpentier, Del Prete
https://arxiv.org/abs/2101.06846

Although the theory is not yet public, my implementation proposes the following improvements:

1. The original formulation only considered linear contact models. If $\delta$ and $\dot{\delta}$ are the penetration depth and the penetration rate along the terrain's normal direction, the new formulation enables using any non-linear contact model belonging to the family of Hunt/Crossley models:

2. The new formulation allows to use non-linear contact models that intrinsically support uneven terrain, under the assumption that we can gather the terrain normal and the terrain height at any $x-y$ coordinate.
3. The new formulation does not integrate the contact dynamics in terms of displacement w.r.t. the anchor point $\mathbf{x} = [\Delta \mathbf{p}; \, \Delta \dot{\mathbf{p}}],$ but it integrates the contact dynamics ignoring the anchor point entirely, therefore $\mathbf{x} = [\mathbf{p}; \, \dot{\mathbf{p}}].$ This enables to do not keep track of the location where the contact was made.
4. The new formulation changes completely how the frictional component of the force (i.e., the tangential component) is computed. It moves the logic to compute the sticking tangential force inside the non-linear contact model, exploiting the fact that we can adopt AD to linearize the contact model with arbitrarily complex logic. The tangential displacement to obtain the $\Delta$ variables are computed by considering an additional state $\mathbf{m}$ of the terrain deformation (that, in a sense, replaces the need of keeping track of the contact creation). This workaround was inspired by the soft-contact model already implemented in the soft contacts of JaxSim. For this reason, the state of the contact dynamics becomes $\mathbf{x} = [\mathbf{p}; \, \dot{\mathbf{p}}; \, \mathbf{m}].$
5. The introduction of $\mathbf{m}$ in the state of the contact dynamics allows us to integrate this variable in continuous time, obtaining a stable and precise computation of the contact forces.
6. The new formulation, thanks to the presence of $\mathbf{m}$ in the state of the contact dynamics and the possibility to have arbitrarily complex non-linear contact models, completely changes the original implementation of sticking/slipping transition. I adapted the sticking/slipping logic used by the default soft contact model of JaxSim to this new contact model, and I obtained sound theoretical results (that also works well in practice, paying the price of a larger $A$ matrix used in $\exp^{At}$).

Open points left to future work:

• The main bottleneck of the whole contact model is the computation of the matrix exponential. In the paper, the authors describe a custom implementation that might be considerably faster than jax.scipy.linalg.expm. However, beyond implementing the logic, which is doable, we should also figure out how to expose AD. It's not clear to me what could be the performance degradation of relying to pure AD w.r.t. the usage of the custom adjoints as done internally in JAX.
• Exposing to the user the logic to enable/disable collidable points is not included in this PR. For the time being, altering directly the low-level boolean array is necessary for thinnin out the considered points.

📚 Documentation preview 📚: https://jaxsim--248.org.readthedocs.build//248/